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Abstract 

This  research  resulted  in  the  development  of  the  Discrete  Element  Method  (DEM) 
as  a  general,  robust,  and  scalable  computer  technique  for  unified  modeling  of  the 
mechanical  behavior  of  solid  and  particulate  materials,  including  the  transition  from  solid 
phase  to  particulate  phase.  Applications  include  gross  damage  of  structures  due  to 
extreme  load  events,  and  high  speed  penetration  of  structures,  involving  materials  such  as 
concrete,  rock,  and  novel  combinations  of  these.  This  computer  approach  is  useful  for 
assessing  vulnerability  of  military  and  civilian  facilities  such  as  nuclear  power  reactors, 
transportation  facilities,  buildings,  and  so  on,  and  for  assessing  efficacy  of  military 
operations  involving  such  structures.  One  of  the  primary  accomplishments  of  this 
research  is  the  development  of  interelement  potentials  for  the  DEM  method  such  that 
accurate  and  convergent  results  are  obtained. 


1.  Executive  Summary 

Although  empirical  equations  for  local  damage  give  a  reasonable  prediction  of 
quantities  such  as  depth  of  penetration,  they  do  not  describe  the  behavior  of  the  concrete 
structure.  In  order  to  improve  the  understanding  of  concrete  subjected  to  severe  loading,  a 
combination  of  experiments  and  numerical  methods  are  needed  for  a  complete  analysis  of 
the  structural  behavior  where  the  material  fracture  mechanisms  can  be  studied  in  detail. 

In  this  work,  a  discrete  element  method  was  developed  as  a  general  computer  technique 
for  unified  modeling  of  the  mechanical  behavior  of  solid  and  particulate  materials, 
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including  the  transition  from  solid  phase  to  particulate  phase.  As  a  result  of  this  research, 
the  following  contributions  in  the  field  of  computer  modeling  by  the  discrete  element 
method  were  made: 


•  A  theoretical  method  was  developed  to  determine  DEM  contact  normal  and  tangential 
stiffnesses  as  functions  of  known  material  properties  so  that  convergence  is  obtained  and 
the  correct  elastic  behavior  is  produced,  depending  on  the  arrangement  of  the  DEM 
elements  in  the  discretization.  For  a  simple  cubic  assembly,  the  normal  and  tangential 
stiffnesses  are 


Kn  =  K,  -  ETt 


(1.1) 


For  a  close-packed  assembly  and  for  an  irregular  assembly,  the  normal  and  tangential 
stiffnesses  are 


K. 
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(1.2) 

(1.3) 


where  E *  and  v*  arc  defined  for  plane  stress  and  plane  strain  situations  in  terms  of  elastic 
modulus  E,  Poisson's  ratio  v,  and  thickness  t  by 
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•  Regarding  contact  failure  for  close-packed  and  irregular  DEM  assemblies,  a  novel 
DEM  contact  behavior  was  developed  in  order  to  account  for  different  failure  modes. 
That  is,  when  failure  is  dictated  by  the  material’s  fracture  toughness  K,c  or  the  material’s 
ultimate  tensile  strength  oull .  The  contact  stress-displacement  curve,  shown  in  Figure 
1.1,  is  proposed  to  be  linear  until  the  contact  normal  and  shear  stresses  reach  critical 
values  given  by 


r  crit 

~  crit  Jn 

”  ~  2  Rt 
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Figure  1.1  DEM  inter-element  constitutive  behavior. 


After  the  contact  forces  reach  the  values  given  by  Equations  1.6  and  1.7,  the  contact 
behavior  is  proposed  to  be  plastic  until  the  area  under  the  stress-displacement  curve  for 
the  normal  direction  spring  reaches  the  material’s  fracture  energy,  given  by 


The  maximum  normal  displacement  is  then  given  by 


8" 

max 


where 
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(1.8) 


(1-9) 


•  For  purposes  of  selecting  a  time  step  size  for  integrating  the  equations  of  motion,  a  new 
method  was  developed  to  accurately  bound  the  maximum  frequency  of  vibration  in  a 
DEM  model,  which  gives  a  more  accurate  prediction  than  the  ad  hoc  approximation  that 
is  usually  implemented. 

•  The  proposed  DEM  method  was  used  to  model  penetration  events  on  concrete  targets. 
To  increase  the  accuracy  of  numerical  analyses  for  projectile  impact,  nonlinear 
compression  behavior  and  strain  rate  effects  were  incorporated  into  the  model. 
Penetration  depth,  scabbing  velocity,  and  ballistic  limit  values  were  obtained  and 
compared  with  available  empirical  formulae  obtained  from  the  literature  and  a  reasonable 
agreement  was  found. 

The  major  findings  are  summarized  below. 

Convergence  Studies 

The  proposed  method  used  a  DEM  model  of  a  unit  cell  of  homogeneous  material 
to  theoretically  develop  criteria  for  contact  stiffness  and  failure  so  that  convergence  of 
elastic  response  and  damage  response  due  crack  growth  is  established  for  the 
representative  models  considered  in  the  study.  The  simulations  showed  very  good 
correlation  (for  both  the  elastic  behavior  and  fracture)  with  theoretical  results.  However, 
the  convergence  rate  was  not  clear  from  the  simulations  performed,  and  it  appears  to 
change  depending  on  element  arrangement  (i.e.  simple  cubic,  close-packed,  or  irregular). 


3 


Nevertheless,  the  use  of  theory  to  determine  contact  parameters  was  very  beneficial 
because  it  avoided  the  tedious  procedure  of  model  calibration  for  a  given  mesh 
refinement. 

Megaclustering  and  Fracture 

The  technique  used  for  solid  model  development  (megaclustering)  proved  to  be 
an  effective  method  of  introducing  irregularity  in  the  model,  with  no  preferential 
directions  for  fracture.  The  technique  allowed  for  the  mesh  to  have  areas  of  stress 
concentration  where  crack  propagation  would  be  more  likely  to  occur.  For  the  type  of 
element  arrangement  produced  by  megaclustering,  the  energy-based  criterion  (influenced 
by  strain  rate  effects)  developed  for  contact  failure  proved  to  be  more  appropriate  than  a 
criterion  based  purely  in  the  material’s  ultimate  tensile  strength  or  fracture  toughness. 

Penetration  Simulations 

The  numerical  results  obtained  in  this  work  showed  that  DEM  is  a  very  useful 
tool  for  simulating  impact  problems  and  projectile  penetration  into  concrete  structures. 
Due  to  the  strain  rates  produced  by  such  events,  material  properties  such  as  tensile 
strength,  compressive  strength,  and  fracture  energy  increased  considerably.  Results  for 
penetration  depth,  scabbing  velocity,  and  ballistic  limit  were  in  agreement  with  the 
bounds  given  by  the  available  empirical  formulae,  especially  the  BRL  equations  to 
predict  scabbing  velocity  and  ballistic  limit.  However,  due  to  the  parameters  chosen  for 
the  model,  global  deformation  and  flexural  cracking  occurred  in  the  models,  which 
influenced  the  level  of  damage  in  the  target. 

Contact  Compression  Behavior 

The  incorporation  of  nonlinear  compression  behavior  to  account  for  material 
“densification”  due  to  concrete  compaction  proved  to  have  a  considerable  influence  on 
the  structural  response  of  a  concrete  target  subjected  to  a  projectile  penetration.  The 
addition  of  this  behavior  allowed  the  medium  to  absorb  energy  in  the  material 
immediately  ahead  of  the  projectile,  allowing  the  numerical  results  to  correlate  better 
with  the  empirical  formulae.  However,  the  implementation  of  this  behavior  increased  the 
ductility  of  the  medium,  influencing  the  amount  of  spalling  at  the  front  face  and  scabbing 
at  the  back  face  of  the  target. 
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Detailed  Summary  of  Work 


2.  Overview  of  the  Discrete  Element  Method  (DEM) 

The  discrete  element  method  (DEM)  discretizes  a  material  using  rigid  elements  of 
simple  shape  that  interact  with  neighboring  elements  according  to  interaction  laws  that 
are  applied  at  points  of  contact.  The  rigid  elements  are  typically  of  circular  or  spherical 
shape  due  to  the  simplicity  and  speed  of  the  contact  detection  algorithm.  However, 
several  researchers  have  applied  the  method  using  elements  of  polygonal  shape  (Heuze  et 
al.  1993;  Kun  and  Herrmann  1996;  Bolander  and  Saito  1997;  Camborde  et  al.  2000; 
Prochazka  2004). 

The  method  is  ideal  for  modeling  many  problems  not  accessible  to  traditional 
continuum-based  methods  such  as  finite  differences  and  finite  elements  since  it  has  the 
advantage  of  inherently  capturing  large  deformations  and  solid  phase  changes  which  are 
commonly  found  in  granular  media.  Other  applications  involving  material  separation  and 
progressive  failure  phenomena  include  concrete  structural  failure,  rock  blasting 
operations,  and  fracture  of  ceramic  or  other  quasibrittle  materials  under  high  velocity 
impact.  In  contrast  to  finite  elements,  the  method  is  particularly  attractive  for  modeling 
geomaterials  due  to  its  ability  construct  a  mesh  that  is  not  completely  continuous  and 
homogeneous,  as  shown  in  Figure  2.1.  Since  the  mesh  is  constructed  by  rigid  elements 
that  interact  with  each  other  at  points  of  contact,  the  DEM  mesh  has  the  capability  of 
easily  constructing  a  medium  with  voids,  imperfections  and  heterogeneities,  as  shown  in 
Figure  2.1(b),  which  are  present  in  rock,  concrete,  asphalt,  and  other  geomaterials. 

a  a 


Figure  2. 1  Finite  element  mesh  (a)  vs.  discrete  element  mesh  (b). 

The  analysis  algorithm  consists  of  three  major  computational  steps:  contact  processing, 
in  which  contact  forces  are  calculated,  element  motion,  in  which  element  displacements 
are  computed,  and  contact  detection,  in  which  new  contacts  are  identified  and  broken 
contacts  are  removed.  In  a  DEM  analysis,  the  interaction  of  the  elements  is  treated  as  a 
dynamic  process  that  alternates  between  the  application  of  Newton’s  second  law  and  the 
evaluation  of  a  force-displacement  law  at  the  contacts.  Newton’s  second  law  gives  the 
acceleration  of  an  element  resulting  from  the  forces  acting  on  it,  including  gravitational 
forces,  external  forces  prescribed  by  boundary  conditions,  and  internal  forces  developed 
at  interparticle  contacts.  The  acceleration  is  then  integrated  to  obtain  the  velocity  and 
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displacement.  The  force-displacement  law  is  used  to  find  contact  forces  from  known 
displacements.  This  process  is  described  in  the  following  sections. 

Element  Interaction  Models 

When  modeling  a  solid  with  no  initial  cracks  (i.e.,  a  nominally  homogeneous 
material),  all  contacts  in  the  assembled  medium  are  initially  bonded  with  the  behavior 
that  is  schematically  illustrated  in  Figure  2.2.  During  the  course  of  a  simulation,  if  a 
bonded  contact  fails,  according  to  criteria  to  be  described,  the  contact  becomes  frictional 
as  shown  also  in  Figure  2.2,  if  indeed  the  contact  still  exists  (i.e.,  if  the  two  DEM 
elements  are  still  being  pressed  into  one  another),  or  the  contact  is  destroyed  if  the  two 
elements  are  separating.  The  computer  implementation  uses  Rayleigh  or  proportional 
damping  (Cook  et  al.  1989),  and  therefore  damping  parameters  are  chosen  as  fractions  of 
critical  damping  at  two  desired  frequencies. 


bonded  contact  frictional  contact 


Figure  2.2  Discrete  element  interaction  models  for  bonded  and  frictional  contacts. 

Bonded  Contact 

In  the  bonded  contact,  the  normal  spring  acts  in  tension  and  compression,  having 
a  linear  elastic  behavior  between  force  and  displacement,  as  shown  in  Figure  2.3(a).  The 
shear  spring  also  has  a  linear  elastic  behavior,  shown  in  Figure  2.3(b),  with  an  elastic 
constant  that  is  usually  lower  than  that  of  the  normal  spring. 


Figure  2.3  Constitutive  relation  for  bonded  contact. 

Frictional  Contact 

The  normal-direction  spring  in  the  frictional  contact  shown  in  Figure  2.4(a)  acts 
only  in  compression.  In  the  tangential  direction,  shown  in  Figure  2.4(b),  the  shear  force 
is  monitored  against  a  maximum  possible  value  according  to  a  Coulomb-type  friction 
law.  The  maximum  shear  force  is  defined  as 
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|/7’/max  [  =  -fx  Fn  +c  (2.1) 

where  fx  and  c  are  the  friction  coefficient  and  the  cohesion  between  the  elements, 
respectively.  If  the  computed  shear  force  F,  is  larger  than  F,Tmx,  F,  is  set  to  F,mm,  allowing 
the  elements  to  slide. 


Figure  2.4  Constitutive  relation  for  frictional  contact. 


Equations  of  Motion 

After  computing  all  the  forces  acting  on  an  element,  the  displacement,  velocities 
and  accelerations  of  each  element  are  computed  using  the  central  difference  method, 
where  details  are  given  in  Tavarez  (2005)  and  Tavarez  and  Plesha  (2006).  Since  explicit 
integration  is  conditionally  stable,  the  time  step  At  must  satisfy 

—  \-l2 -%)  (2.2) 

"max 

for  linear  viscous  damping,  where  §  is  the  fraction  of  critical  damping  at  ft>max,  which  is 
the  highest  natural  frequency  of  the  mesh.  For  practical  problems,  c is  rarely  known 
because  it  is  expensive  to  determine.  Therefore,  an  accurate  upper  bound  for  this 
frequency  is  desirable. 


Frequency  Bound 

In  this  work,  we  have  developed  a  new  approach  for  obtaining  an  upper  bound  for 
cumax  by  computing  the  largest  eigenvalue  for  a  unit  cell  of  DEM  elements  as  shown  in 
Figure  2.5  This  unit  cell  consists  of  seven  rigid  elements  arranged  in  a  close-packed 
assembly.  A  complete  mesh  for  a  structure  is  then  constructed  by  discretizing  the 
medium  using  a  large  number  of  these  DEM  unit  cells  inter-connected  with  one  another. 
Therefore,  (&>max  )e  can  now  be  obtained  by  calculating  the  largest  eigenvalue  of  the  DEM 
unit  cell  structure.  As  will  be  shown  in  Section  3,  the  construction  and  analysis  of  this 
DEM  unit  cell  was  also  necessary  to  obtain  the  correct  values  for  contact  stiffness  for 
modeling  elastic  material  behavior. 

Mass  and  stiffness  matrices  were  computed  for  the  DEM  unit  cell  shown  in 
Figure  2.5,  which  consists  of  21  degrees  of  freedom  (14  translational  d.o.f.  and  7 
rotational  d.o.f.).  Values  for  the  stiffness  and  mass  matrices  are  given  in  Appendix  A  of 
Tavarez  (2005).  Frequency  analysis  of  the  unit  cell  shown  in  Figure  2.5  provides 

("maxi  @  3.742^/A^  jm  ,  where  m  is  the  mass  of  each  DEM  element  and  Kn  is  the  normal 
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contact  stiffness.  Substituting  this  frequency  into  Equation  2.2,  we  obtain  an  accurate 
upper  bound  estimate  for  the  critical  time  step  A tcrU  as 


A',„ -0.535^VX  (/ThF-^ 


(2.3) 


Figure  2.5.  DEM  discretization  and  unit  cell  used  for  computing  maximum  element 
frequency. 


3.  DEM  for  modeling  brittle  materials 

We  use  the  concept  of  clusters  to  represent  solid  media.  Clusters  are  formed  by 
bonding  a  number  of  circular-shaped  DEM  elements  into  a  semi-rigid  configuration.  To 
model  a  solid  such  as  concrete,  this  concept  can  be  extended  to  megaclustering  wherein 
the  entire  volume  of  a  solid  is  modeled  by  bonding  individual  clusters  together.  Once  the 
model's  geometry  is  established,  boundary  conditions  are  prescribed.  This  generation 
scheme  will  allow  for  "molding"  of  structures  or  structural  components  of  arbitrary 
shape. 

Determination  of  Inter-Element  Normal  and  Shear  Stiffness 

The  macroscopic  linear  elastic  behavior  of  isotropic  materials  can  be 
characterized  by  two  elastic  constants:  Young’s  modulus  E  and  Poisson’s  ratio  v,  and  the 
physical  parameters  governing  DEM  element  interactions  have  typically  been  determined 
by  the  ad  hoc  process  of  validating  a  numerical  simulation  of  a  standard  laboratory  test 
with  an  actual  experimental  result  (Potyondy  et  al.  1996;  Magnier  and  Donze  1998). 
Unfortunately,  if  a  new  discretization  for  a  particular  problem  is  produced  but  with 
different  DEM  element  size  (such  as  making  element  size  smaller  to  improve  accuracy  of 
the  simulation),  the  calibration  process  would  need  to  be  repeated  in  its  entirety  to  obtain 
the  new  DEM  parameters. 

In  our  study  we  developed  a  new  approach  wherein  we  theoretically  established 
the  interelement  normal  and  tangential  stiffnesses  (and  failure  parameters,  to  be  discussed 
later)  as  functions  of  element  size  and  commonly  accepted  elastic  constants  including 
Young's  modulus  and  Poisson's  ratio.  This  was  done  by  the  analysis  of  a  DEM  unit  cell  of 
material.  Figure  3.1(a)  shows  an  isotropic  solid  material  element  (with  known  elastic 
modulus  and  Poisson’s  ratio)  subjected  to  uniaxial  stress.  The  volume  of  material  is  then 
modeled  using  the  DEM  unit  cell  shown  in  Figure  3.1(b).  Using  this  unit  cell,  the  normal 
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and  tangential  spring  constants  can  then  be  determined  as  functions  of  the  known  elastic 
modulus  and  Poisson’s  ratio  of  the  material,  and  if  needed  the  DEM  element  size. 


(«) 
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2  R  sin  — 
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2  R 


2  R  sin  —  - 
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tttttttttt 


R 
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Figure  3.1  DEM  unit  cell  for  determination  of  inter-element  spring  constants. 


The  unit  cell  in  Figure  3.1(b)  contains  7  elements  having  3  degrees  of  freedom  per 
element  (two  translations  and  one  rotation).  By  requiring  the  deformation  of  this  unit  cell 
to  agree  with  the  exact  response  for  an  isotropic  material  (see  Tavarez  2005),  we  obtain 
the  normal  and  tangential  stiffnesses  to  be 


Kn  =  -r =r - \E't 

■  V3(l-v-) 


K, 


.iziA . 

(1+v  )  yfl\-V  ) 


(3.1) 

(3.2) 


where  E *  and  v*  are  defined  for  plane  stress  and  plane  strain  situations  by 


plane  stress 

plane  strain 

E 

E 

7i — n 

(3.2) 

(l-v  ) 

V 

V 

(l-v) 

(3.4) 

DEM  Fracture  Criteria 

Different  approaches  have  been  used  in  the  past  to  model  fracture  in  DEM,  where 
most  of  this  work  considers  a  tensile  force  limit  for  contact  failure  (Zubelewicz  and 
Bazant  1987;  Masuya  et  al.  1994;  Bolander  and  Saito  1997;  Brara  et  al.  2001;  Mishra  and 
Thorton  2001).  That  is,  if  the  tensile  force  in  a  certain  contact  exceeds  a  limit  value,  the 
contact  breaks  and  can  no  longer  support  tensile  forces.  This  force  limit  is  usually 
determined  by  calibrating  DEM  simulation  results  to  available  experimental  data  from 
tensile  and  compressive  tests  performed  on  a  particular  material.  However,  since  force 
transmission  through  a  DEM  medium  can  only  occur  via  the  inter-element  contacts,  the 
number  of  contacts  within  a  DEM  medium  will  most  likely  affect  the  fracture  behavior. 
Therefore,  this  contact  force  limit  is  most  likely  to  be  element  size  dependent.  Taking 
into  account  this  element  size  dependency,  several  researchers  have  considered  a  failure 
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criterion  based  on  the  tensile  strength  of  the  material,  using  a  DEM  model  consisting  of 
two  elements,  as  shown  in  Figure  3.2  (Sawamoto  et  al.  1998).  The  maximum  normal 
stress  failure  criterion  states  that  failure  occurs  when  the  maximum  principal  stress 
reaches  the  tensile  strength  of  the  material.  The  tensile  force  required  to  break  the 
contact  is  then  determined  by 

(3.5) 

where  aull  is  the  tensile  strength  of  the  material.  As  seen  in  Equation  3.5,  this  criterion 
suggests  that  the  tensile  force  required  for  contact  failure  is  linearly  proportional  to  the 
DEM  element  radius  R. 


Figure  3.2  Tensile  strength  criterion  for  fracture. 

Several  researchers  have  imposed  a  limit  in  both  the  tensile  and  shear  stress  in  the 
contact,  as  shown  in  Figure  3.3  (Potyondy  et  al.  1996;  Donze  et  al.  1997;  Magnier  and 
Donze  1998).  These  thresholds  are  also  selected  according  to  shear  and  tensile  strengths 
determined  in  laboratory  experiments. 


/, 


j'  max 

j'max 

f  compression 

J  n 

*  tension 

^  J n  ^  | 

Figure  3.3  Failure  surface  for  bonded  contact  considering  tensile  and  shear  limits. 
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Fracture  Toughness  Failure  Criterion 

In  what  follows,  we  theoretically  develop  a  contact  failure  criterion  based  on 
DEM  element  size  and  known  material  properties.  Since  this  method  will  eventually  be 
used  to  model  extensive  cracking  and  fragmentation  in  a  medium,  a  criterion  based  only 
on  the  ultimate  tensile  strength  of  the  material  would  produce  inaccurate  results  when 
stress  singularities  arise  due  to  flaws  and  cracking  in  the  medium.  The  following 
approach  addresses  this  issue. 

The  results  to  be  described  here  begin  with  the  original  work  of  Potyondy  and 
Cundall  (2004).  In  their  work,  they  postulate  a  breaking  strength  between  the  bonds  of 
contacting  DEM  elements,  and  then  address  the  issue  of  determining  the  fracture 
toughness  of  the  material  that  is  being  modeled  by  the  DEM  discretization.  The  approach 
described  herein  differs  from  theirs  in  that  we  view  the  fracture  toughness  as  a 
fundamental  material  parameter,  and  then  attempt  to  use  theory  to  develop  expressions 
for  the  behavior  of  bonds  between  contacting  DEM  elements.  A  main  feature  of  such  an 
approach  is  that  it  promises  to  provide  convergence.  That  is,  it  provides  a  precise 
specification  of  how  the  contact  behavior  between  DEM  elements  must  change  as  a 
function  of  element  size  so  that  convergence  to  the  exact  solution  is  achieved  as  the  DEM 
element  size  becomes  smaller  during  model  refinement.  To  model  the  progressive  failure 
of  a  solid  due  to  crack  growth,  a  local  rupture  criterion  is  applied  to  the  bonded  contacts 
between  interacting  elements.  The  plane  stress/strain  near  crack  tip  stress  field  for  pure 
Mode  I  loading  in  a  homogeneous,  isotropic,  linear  elastic  material  is  given  by  (Anderson 
1994): 


y/2n r 


A„r"nF<"'(e) 


(3.6) 


Since  crack  growth  is  controlled  by  the  stresses  near  the  crack  tip,  i.e.,  at  small  r,  the  first 
term  in  Equation  3.6  gives  the  highest  contribution  to  the  stress  field  at  the  crack  tip. 
Therefore,  keeping  only  the  first  term  in  Equation  3.6,  the  well  known  expression  for  the 
stresses  close  to  the  crack  tip,  shown  in  Figure  3.4,  is  obtained. 


Figure  3.4  Crack  tip  stress  field  for  pure  Mode  I  cracking. 


Restricting  attention  to  8  =  0,  the  stress  aee  is  then 


K, 

V 2nr 


(3.7) 
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The  force  required  to  advance  the  crack  by  breaking  the  contact  at  the  crack  tip  of  a  DEM 
model  is  determined  by  integrating  the  stress  field  over  the  DEM  element  diameter  and 
thickness  (Potyondy  and  Cundall  2004),  as  shown  in  Figure  3.4.  Therefore, 


/  2  R  t  2  R 

f  =  J  J  amdrdt  =J  j  ’tr~Kndrdt 
oo  oo  V  2/r 


(3.8) 


(3.9) 


V  n 

where /is  defined  in  Figure  3.2.  According  to  Irwin’s  criterion,  cracking  will  occur  when 
the  value  of  K,  reaches  the  material's  fracture  toughness  Kk.  Therefore,  the  required  force 
to  break  the  contact  is  given  by 


2kJ~  <3A0> 

V  K 

As  shown  in  Equation  3.10,  material  failure  parameters  in  this  approach  are  dependent  on 
DEM  element  size  and  the  fracture  toughness  of  the  material.  As  described  above,  this 
result  is  very  significant  because  it  prescribes  exactly  how  the  DEM  interelement  strength 
must  change  as  a  function  of  element  size  so  that  convergence  is  obtained  as  element 
radius  R  vanishes  in  the  limit  of  mesh  refinement. 


4.  Convergence  of  DEM  simulations 

Before  discussing  convergence  of  DEM  simulations,  it  is  instructive  to  briefly 
discuss  convergence  of  simulations  performed  using  the  finite  element  method  (FEM). 
In  a  typical  FEM  simulation,  a  problem  is  discretized  using  a  user-selected  number  of 
degrees  of  freedom.  If  the  problem  being  modeled  is  elliptic  (linear-elastic,  small 
deformations,  no  singularities  such  as  point  forces  or  cracks,  etc.)  then  it  is  expected  that 
accuracy  of  results  will  improve  as  the  model  is  refined  by  adding  more  degrees  of 
freedom,  and  that  convergence  to  the  exact  solution  will  be  obtained  in  the  limit  that  an 
infinite  number  of  degrees  of  freedom  are  used.  Fortunately  and  very  importantly,  there 
are  theoretical  results  showing  that  FEM  simulations  indeed  do  converge,  and 
furthermore  the  rate  of  convergence  for  most  popular  element  types  is  known.  In  this 
study  we  attempt  to  help  establish  similar  results  for  DEM  simulations. 

We  have  studied,  and  established  the  convergence  properties  for  DEM  models  for 
applications  to  axial  deformation  of  a  bar,  including  the  Poisson  effect,  and  bending 
deformation  of  a  beam,  for  simple  cubic  arrangements  of  DEM  elements,  close  packed 
arrangements,  and  irregular  (random)  arrangements.  In  all  cases,  the  interelement 
potentials  we  have  described  in  Section  3  give  convergence  as  the  number  of  DEM 
element  used  in  the  discretization  is  increased.  In  this  section,  we  present  results  for  just 
one  of  these  studies,  namely  bending  deformations  using  a  simple  cubic  arrangement  of 
DEM  elements. 

Convergence  for  Elastic  Response  of  a  Cantilever  Beam 

A  DEM  model  of  a  cantilever  beam  subjected  to  a  sinusoidal  tip  load  with 
frequency  co  was  considered.  A  closed-form  solution  for  this  problem  can  be  obtained  by 
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idealizing  the  structure  as  a  single  degree  of  freedom  (SDOF)  model  with  natural 
frequency  equal  to  the  fundamental  natural  frequency  of  a  continuous  cantilever  beam 

(i.e.,  (jon  =  3 ,52-yj El / ML7,  ).  However,  this  solution  assumes  that  only  the  fundamental 

mode  shape  of  the  cantilever  beam  is  excited  by  the  loading  conditions.  Since  higher 
modes  might  not  be  insignificant  for  the  degree  of  precision  needed,  the  accuracy  of  the 
DEM  simulations  for  this  problem  will  be  studied  by  also  comparing  the  DEM  results 
with  simulations  performed  using  finite  element  analysis  (FEA).  The  finite  element 
analysis  was  performed  in  ANSYS,  where  the  cantilever  beam  was  discretized  using 
1440  four-node  plane  quadrilateral  elements. 

Figure  4. 1  shows  three  DEM  discretizations  for  the  cantilever  beam  using  close- 
packed  assemblies  with  81,  729,  and  6561  elements,  respectively.  The  element  size  was 
reduced  by  a  factor  of  3  for  each  refinement.  All  DEM  simulations  use  the  values  for 
normal  and  tangential  spring  stiffness  given  by  Equations  3.1  and  3.2  with  E  =  1.2  GPa 
and  v  =  0.3.  The  forcing  frequency  was  chosen  to  be  0.2  times  the  value  of  the 
fundamental  natural  frequency,  and  a  fraction  of  critical  damping  of  §  =  10  %  was  used 
at  this  frequency.  Figure  4.2(a)  shows  the  tip  displacement  as  a  function  of  time  for  the 
three  models.  Figure  4.2(b)  shows  the  tip  displacement  obtained  by  the  finest  mesh  along 
with  the  SDOF  solution  and  the  FEA  solution. 


Figure  4.1  DEM  models  of  cantilever  beam  using  a  close-packed  assembly. 

As  shown  in  Figure  4.2(b),  there  is  also  excellent  agreement  between  the  DEM  solution 
and  both  the  SDOF  solution  and  the  FEA  solution,  and  results  also  appear  to  converge  to 
a  definite  value  in  the  limit  of  mesh  refinement.  However,  since  the  ratio  (u/mn  was 
chosen  to  be  greater  for  these  models  ( (o/a)n  =0.2 ),  the  importance  of  the  higher  modes 
of  vibration,  which  are  neglected  in  the  SDOF  model,  is  greater  and  thus,  the  DEM 
solution  agrees  slightly  better  with  the  finite  element  solution.  Applying  Richardson 
extrapolation  (Cook,  et  al  1989)  to  the  results  of  these  three  DEM  models  provides  an 
estimate  of  8nmx  =  0.4754  m  as  the  steady  state  maximum  amplitude  of  vibration  that  the 
DEM  models  will  converge  to.  The  finite  element  solution  provides 
8 FEM  =  0.4749  m  which  is  within  0.1  %  of  the  DEM  extrapolated  solution. 
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Figure  4.2  DEM  tip  displacement  results  and  comparison  with  theoretical  solution  using 
a  close-packed  assembly. 


We  also  conducted  convergence  analyses  for  irregular  (random)  DEM  discretizations. 
While  these  results  appear  to  show  convergence,  it  is  interesting  to  note  that  contrary  to 
the  cases  using  simple  cubic  and  close-packed  arrays,  the  irregularly  structured  models 
seem  to  become  more  flexible  as  the  mesh  is  refined.  Moreover,  results  appear  to 
converge  to  a  solution  that  is  lower  (about  6.5%)  than  the  FEA  solution.  Therefore,  it  is 
concluded  that  when  Equations  3.1  and  3.2  are  used  for  an  irregular  DEM  assembly, 
results  may  underpredict  the  steady  state  response  of  the  beam.  One  possible  explanation 
for  why  the  models  become  stiffer  as  the  DEM  mesh  is  refined  (for  the  simple  cubic  and 
close-packed  arrangements)  could  be  that  for  each  refinement,  there  are  contact  springs 
that  are  closer  to  the  top  and  bottom  surfaces  of  the  beam,  increasing  the  flexural  stiffness 
of  the  model.  This  appears  to  have  a  greater  effect  on  the  behavior  of  the  model  than 
adding  more  degrees  of  freedom,  which  usually  results  in  making  a  model  more  flexible. 
On  the  other  hand,  the  addition  of  degrees  of  freedom  appears  to  have  a  greater  effect  on 
the  behavior  of  the  model  when  using  an  irregular  arrangement,  resulting  in  a  numerical 
solution  that  becomes  “softer”  as  the  mesh  is  refined. 


Rate  of  Convergence 

The  purpose  of  having  at  least  three  different  results  for  each  element 
arrangement  pattern  was  to  determine  the  rate  of  convergence  for  the  DEM  analysis. 
Using  the  results  of  three  meshes,  the  rate  of  convergence  (if  uniform)  can  be  computed 
using  Richardson  extrapolation  (Cook  et  al.  1989)  as 


P  =  log 


0|  02 
02  —  03  t 


1 

log  a 


(4.1) 


where  are,  in  this  case,  the  values  for  the  maximum  steady  state  displacement  for  the 
three  different  discretizations,  and  a  is  the  element  size  reduction  factor  for  each  mesh 
refinement,  which  must  be  a  constant  for  each  refinement.  Table  4.1  shows  the  rate 
convergence  computed  for  the  three  different  meshes  used  for  representing  the  cantilever 
beam. 
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Table  4. 1  Summary  of  convergence  results  for  the  three  model  lattices  used. 


DEM  Element 
Arrangement 

Number  of  Particles  Used 

a 

m 

Mesh  1 

Mesh  2 

Mesh  3 

simple  cubic 

90 

810 

7290 

3 

4.57 

close-packed 

81 

729 

6561 

3 

1.52 

irregular 

600 

1200 

2400 

ri 

3.00 

Unfortunately,  the  rate  of  convergence  was  different  for  each  element  arrangement 
pattern.  As  shown  in  this  table,  the  simple  cubic  element  arrangement  pattern  showed  the 
largest  rate  of  convergence  with  a  value  of  p  =  4.57  .  However,  this  surprisingly  large 
value  is  no  doubt  not  accurate  and  might  be  due  to  the  fact  that  the  use  of  only  90 
elements  in  a  simple  cubic  discretization  is  too  coarse  of  a  mesh  for  this  problem  to 
provide  a  uniform  rate  of  convergence.  Consequently,  the  result  of  mesh  1  is 
considerably  lower  than  the  results  of  meshes  2  and  3,  yielding  an  apparent  large  rate  of 
convergence.  The  same  argument  can  be  conjectured  for  the  rate  of  convergence 
computed  for  the  irregular  discretization.  On  the  other  hand,  the  rate  of  convergence  for 
the  close-packed  discretization  was  p  =  1.52,  which  is  a  more  reasonable  value. 
However,  while  these  results  are  the  first  of  their  kind  that  we  are  aware  of,  considerable 
additional  work  is  still  required  to  adequately  understand  the  convergence  properties  of 
DEM  simulations  of  elastic  response  in  solids. 

Convergence  for  Damage  Response  due  to  Crack  Growth 

In  order  to  study  convergence  when  the  response  includes  damage,  three  DEM 
models  representing  a  plate  with  a  crack  subjected  to  pure  Mode  I  loading  were 
developed.  The  material  modeled  by  these  DEM  meshes  has  Young’s  modulus 
E  =  29.1  GPa,  Poisson's  ratio  v  =  0.3,  and  fracture  toughness  Klc  =1.0MPaVm  .  The 

models,  shown  in  Figure  4.3,  consist  of  256,  1024,  and  4096  elements,  respectively.  As 
such,  the  DEM  element  radius  was  reduced  by  a  factor  of  2  with  each  mesh  refinement. 
The  edge  crack  (shown  by  the  solid  line  at  the  center  of  each  model),  was  created  by 
giving  contacts  along  the  crack  the  standard  behavior  described  by  the  frictional  contact 
mechanism  shown  in  Figure  2.4.  That  is,  if  two  initially  contacting  DEM  elements  on  the 
crack  surface  are  separated  (which  is  the  case  for  all  of  the  simulations  of  the  models  in 
Figure  4.3),  there  are  no  forces  acting  between  them,  while  if  they  remain  in  contact,  they 
have  Coulomb  friction  response. 

As  seen  in  Figure  4.3,  these  models  were  constructed  using  a  simple  cubic 
assembly  of  elements.  This  type  of  assembly  was  chosen  due  to  its  simplicity  for 
theoretical  analysis  leading  to  Equation  3.10.  Even  though  symmetry  could  have  been 
used  for  this  type  of  problem  to  reduce  the  number  of  elements,  it  was  decided  to  model 
the  entire  problem  in  order  to  avoid  the  application  of  symmetry  constraints  to  each  of  the 
DEM  elements  located  at  the  symmetry  boundary.  For  the  loading  shown  in  Figure  4.3, 
the  displacements  are  such  that  the  crack  surface  undergoes  zero  vertical  displacement, 
while  the  upper  boundary  moves  in  the  upward  direction  and  the  lower  boundary  moves 
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an  equal  amount  in  the  downward  direction.  With  a  simple  cubic  assembly,  the  crack 
plane  is  perpendicular  to  the  direction  of  the  far  field  stress  and  thus  the  force  at  the  crack 
tip  is  carried  by  only  the  normal  direction  springs  ahead  of  the  crack  tip.  For  these 
models,  the  normal  and  tangential  stiffnesses  were  assumed  to  be 
Kn=K,=  Et  (4.2) 

The  results  of  the  DEM  simulations  are  shown  in  Figure  4.3,  where  the  solid  data  points 
represent  values  of  the  displacement  and  far  field  stress  that  cause  the  crack  to  propagate 
unstably  across  the  model.  Very  importantly,  these  results  appear  to  converge  and  an 
extrapolated  value  for  the  far  field  stress  that  causes  unstable  crack  growth  can  be 
computed.  Applying  Richardson  extrapolation  to  the  results  of  these  three  models 
provides  an  estimate  of  ouh  =  2.78  MPa  as  the  value  of  the  far  field  stress  that  the  DEM 
models  will  converge  to. 
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Figure  4.3  Convergence  for  crack  growth  in  a  DEM  model  of  a  cracked  plate  subjected 


to  pure  Mode  I  loading;  solid  circles  in  the  plot  represent  the  value  of  far  field  stress  and 


displacement  that  cause  unstable  crack  growth. 


It  is  important  to  mention  that  an  exact  analytic  solution  for  this  problem  (finite 
width  plate  with  an  edge  crack)  does  not  exist.  This  is  due  to  the  difficulty  of  finding  a 
solution  that  satisfies  the  traction-free  boundary  conditions  on  the  free  edges.  However, 
several  numerical  techniques  have  been  developed  in  order  to  obtain  very  accurate 
approximate  solutions  to  this  type  of  problem.  Such  techniques  include  the  Boundary 
Collocation  Method,  Green’s  Function  Method,  Asymptotic  Approximation,  and  the 
Finite  Element  Method,  among  others.  A  commonly  used  solution  to  this  problem, 
consisting  of  an  finite  width  plate  with  an  edge  crack  subjected  to  pure  Mode  I  loading  is 
given  by  (Gross  and  Strawley,  1964;  Brown  and  Strawley,  1966) 


^.12-0.23(a/c)+10.6(a/c)2  -21.7(a/c)3  +30.4(o/c)’]1 


16 


where  oF  is  the  value  of  the  far  field  stress  that  causes  the  crack  to  propagate  unstably 
across  the  plate,  a  is  the  crack  length  (=  0.03125  m)  and  c  is  the  width  of  the  plate 
(a/c  =  1/8  for  all  the  models  of  Figure  4.3).  Equation  4.3  was  obtained  by  a  least  square 
fitting  to  data  obtained  by  the  Boundary  Collocation  Method,  which  is  accurate  to  within 
0.5  %  for  a/c  <,  0.6 .  Using  the  parameters  cited  earlier,  Equation  4.3  provides 
af  =  2.61  MPa  .  The  extrapolated  DEM  solution  overpredicts  the  solution  in  Equation 
4.3  by  about  6.5  %.  While  part  of  this  disagreement  may  be  due  to  error  in  the 
Richardson  extrapolation  and  the  error  in  Equation  4.3  itself,  it  is  nonetheless  clear  that 
the  DEM  results  are  converging  to  a  value  that  is  greater  than  the  exact  solution.  As 
mentioned  before,  the  value  of  the  tangential  spring  stiffness  influences  the  result,  and  it 
is  likely  that  the  use  of  Equation  4.2  for  the  tangential  stiffness  is  inadequate.  Another 
potential  reason  for  this  disagreement  involves  the  incorrect  fracture  energy  release  rate 
obtained  by  using  the  present  approach  for  contact  stiffness  and  failure  criterion.  This 
will  be  explained  in  more  detail  in  Section  5.  While  we  view  the  agreement  between 
DEM  and  the  exact  solution  as  good,  additional  work  was  needed  in  order  to  improve  the 
accuracy  of  the  DEM  results.  Nevertheless,  the  proposed  failure  criterion  based  on  the 
fracture  toughness  of  the  material  is  superior  and  better  suited  for  problems  with  existing 
cracks  than  a  criterion  based  on  the  material’s  tensile  strength,  the  results  of  which  fail  to 
converge. 


5.  General  criterion  for  DEM  interelement  failures 

In  Section  3,  a  criterion  for  contact  failure  was  developed  as  a  function  of  element 
size  and  the  material’s  fracture  toughness,  and  convergence  results  presented  in  Section  4 
revealed  that  for  a  model  of  a  plate  with  an  edge  crack  (using  a  simple  cubic  assembly), 
this  criterion  produced  results  which,  while  in  very  good  agreement  with  the  exact 
solution,  they  appeared  to  converge  to  a  far  field  stress  that  was  slightly  higher.  This 
section  describes  an  enhanced  failure  criterion  that  appears  to  provide  convergence,  while 
allowing  for  failure  according  to  the  maximum  principal  stress  failure  criterion  for 
uncracked  media,  and  failure  according  to  fracture  mechanics  for  media  with  cracks. 

Maximum  Principal  Stress  Failure  Criterion 

In  a  DEM  model  with  no  initial  flaws,  where  the  fracture  toughness  is  a  constant 
material  property,  a  failure  criterion  based  on  fracture  mechanics  and  the  material’s 
fracture  toughness  will  be  inappropriate  and  will  overpredict  a  structure's  strength 
compared  to  the  maximum  principal  stress  criterion.  Due  to  the  similarity  of  the  discrete 
element  method  with  molecular  dynamics,  the  fracture  toughness  will  vary  depending  on 
the  molecular  (element)  lattice.  As  such,  since  there  are  no  guaranties  that,  in  the  limit  of 
mesh  refinement,  that  a  DEM  model  will  have  flaws  greater  than  the  critical  crack  size 
for  the  material  being  modeled,  hence  we  have  developed  an  alternate  criterion  that 
automatically  allows  failure  to  be  dictated  by  the  ultimate  tensile  strength  of  the  material, 
ouU ,  or  the  material’s  fracture  toughness,  Klc,  as  appropriate. 
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Fracture  Energy  Failure  Criterion 

Figure  5.1  shows  a  cracked  specimen  before  and  after  an  infinitesimal  crack 
extension  A  a.  The  energy  release  rate  in  a  homogeneous,  isotropic  linear  elastic  material 
is  determined  by  calculating  the  total  work  balance  per  unit  thickness  for  extending  a 
crack  a  distance  A  a.  Therefore, 
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where  df  is  the  crack  opening  displacement  at  position  x  behind  the  crack  tip  after  crack 
extension  Aa,  and  oyy  is  a  function  of  both  position  x  and  the  crack  opening  displacement 
6  behind  the  crack  tip. 
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Figure  5.1  Specimen  configurations  before  and  after  crack  advance  A  a. 


Solving  Equation  5.1  and  taking  the  limit  as  Aa  -»  0  yields  the  well  known  relationship 
between  energy  release  rate  and  stress  intensity  factor  K, 


In  the  present  work,  it  will  be  assumed  that  the  energy  required  to  fracture  an 
isotropic  material  is  independent  of  the  mode  of  fracture.  Therefore,  the  fracture  energy 
Gf  or  critical  energy  release  rate  is  then  given  by  Equation  5.2  when  the  value  of  K, 
reaches  the  material's  fracture  toughness  Kk.  That  is, 


(5.3) 


In  DEM,  Equation  5.1  simplifies  since  A  a  is  equal  to  the  element  diameter  and 
therefore  oyy  is  averaged  by  dividing  the  contact  force  by  the  diameter  and  thickness  of 
the  DEM  element.  Moreover,  the  crack  opening  displacement  is  also  averaged  over  the 
diameter  of  the  element  and  consists  of  the  stretch  of  the  normal  direction  spring  in  the 
crack  tip  contact.  Therefore,  in  the  present  approach,  the  critical  energy  release  rate  is 
obtained  by  computing  the  area  under  the  curve  of  ayy  vs.  <5  at  contact  failure.  For 
instance,  Figure  5.2  shows  the  oyy  vs.  6  behavior  when  Equation  3.10  is  used  as  a 
criterion  for  normal  contact  failure,  and  computes  the  fracture  energy  for  this  condition. 
As  shown  in  the  analysis  in  Figure  5.2,  for  a  simple  cubic  assembly  ( Kn  =Et),  the 


obtained  fracture  energy  will  be 


Gf  = 


Kl 
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which  is  different  from  the  result  in  Equation  5.3.  More  importantly,  this  result  shows 
that  in  order  to  have  the  correct  critical  energy  release  rate  using  Equation  3.10,  the 
normal  contact  stiffness  at  the  crack  tip  must  be  reduced  by  a  factor  of  rc. 


Figure  5.2  DEM  Fracture  energy  computation  using  Equation  3.10  for  contact  failure. 

Another  important  issue  to  consider  is  that  the  criterion  in  Equation  3.10  considers  only 
Mode  I  cracking,  where  there  is  no  force  being  carried  by  the  tangential  spring  at  the 
crack  tip  contact.  Since  in  a  two-dimensional  simulation  both  Modes  I  and  II  can  be 
present  simultaneously  in  a  contact,  especially  for  close-packed  and  irregular 
discretizations,  a  fracture  criterion  should  be  developed  to  account  for  this  type  of  loading 
where  both  the  normal  and  tangential  springs  at  the  crack  tip  contact  are  being  deformed. 

In  what  follows,  both  a  DEM  flaw-free  medium  and  a  DEM  cracked  medium  will 
be  considered  in  order  to  develop  a  combined  criterion  for  contact  failure  where  the 
ultimate  tensile  strength  of  the  material  is  used  together  with  the  material’s  fracture 
toughness  to  model  crack  tip  contact  failures. 


DEM  Flaw-Free  Medium 

Returning  to  the  DEM  unit  cell,  shown  in  Figure  5.3,  we  can  obtain  both  the 
normal  and  tangential  contact  forces  as  a  function  of  the  applied  uniaxial  tensile  stress  in 
the  unit  cell.  When  the  applied  stress  approaches  the  material’s  tensile  strength  oul„  using 
the  unit  cell  results  for  normal  and  tangential  contact  stiffness,  these  forces  are 


(5.5) 


(5.6) 


These  critical  forces  for  a  close-packed  assembly  are  analogous  to  the  criterion  leading  to 
Equation  3.5,  which  is  only  based  on  the  tensile  strength  of  the  material. 

For  a  flaw-free  medium,  contact  normal  and  tangential  forces  will  reach  the 
critical  forces  given  by  Equations  5.5  and  5.6  simultaneously  and  then  the  model  will 
reach  the  ultimate  tensile  strength  of  the  material.  Regardless  of  the  material’s  fracture 
toughness,  the  contact  forces  should  not  exceed  those  given  by  Equations  5.5  and  5.6. 
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Otherwise,  the  model’s  tensile  strength  will  be  greater  than  the  material’s  tensile  strength. 
However,  a  DEM  cracked  medium  presents  a  different  scenario  for  ultimate  failure,  as 
follows. 


Figure  5.3  DEM  unit  cell  used  for  determination  of  critical  contact  forces. 

DEM  Cracked  Medium 

When  cracks  are  present,  the  contact  forces  at  the  crack  tip  reach  the  critical 
forces  given  by  Equations  5.5  and  5.6  relatively  quickly  during  loading  since  stresses  at  a 
DEM  crack  tip  become  infinite  in  the  limit  of  mesh  refinement.  However,  it  is  proposed 
that  complete  failure  should  not  occur  before  the  contact  achieves  the  material’s  fracture 
energy  Gf .  In  order  to  satisfy  this  condition,  the  stress-displacement  curve  shown  in 
Figure  6.7  is  adopted  for  both  the  normal  and  tangential  contact  at  the  crack  tip. 


Figure  5.4  DEM  contact  stress-displacement  behavior  and  fracture  energy. 


As  shown  in  Figure  5.4,  the  contact  springs  (both  normal  and  tangential)  behave  linearly 
up  to  the  point  where  the  forces  satisfy  Equation  5.5  or  Equation  5.6  for  the  normal  and 
tangential  springs,  respectively.  Once  these  forces  are  reached,  the  behavior  is  proposed 
to  be  plastic  until  the  area  under  the  stress-displacement  curve  for  the  normal  spring, 
shown  in  Figure  5.4(a),  reaches  the  fracture  energy  of  the  material,  given  by  Equation 
5.3.  It  should  be  noted  that  the  area  under  the  curve  in  Figure  5.4(a)  is  the  minimum 
energy  lost,  as  energy  may  also  be  stored  in  the  tangential  direction  of  the  contact. 
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The  maximum  normal  displacement  <5”ax  is  obtained  by  computing  the  area  under 
the  curve  in  Figure  5.4(a)  and  equating  this  to  the  material’s  fracture  energy  Gf. 
Therefore,  the  maximum  normal  displacement  is  given  by 
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Therefore,  Equation  5.7  gives  a  displacement-based  failure  criterion,  where  the  maximum 
displacement  is  a  function  of  the  material’s  tensile  strength,  elastic  modulus,  Poisson’s 
ratio,  and  fracture  toughness.  Since  contact  stresses  are  computed  using  the  radius  of  the 
elements  in  contact,  the  failure  criterion  is  also  element  size  dependent.  It  should  be 
noted  that  if  unloading  occurs  in  the  contact  after  it  has  reached  the  plastic  range,  the 
slope  of  the  unloading  force-displacement  curve  will  be  equal  to  the  initial  slope  ( Kn  and 
K,  for  the  normal  and  tangential  stiffnesses,  respectively). 


Convergence 

To  test  the  failure  criterion  given  in  Equation  5.7,  three  DEM  models  representing 
a  plate  with  a  crack  subjected  to  pure  Mode  I  loading  were  developed.  The  models, 
shown  in  Figure  5.5,  consist  of  1 172,  4713,  and  18899  elements,  respectively,  arranged  in 
a  close-packed  assembly.  The  material  modeled  by  these  DEM  meshes  has  elastic 
modulus  E  =  29.7  GPa,  Poisson's  ratio  v  =  0.3,  ultimate  tensile  strength  ouh  =  3.44  MPa , 

and  fracture  toughness  KIc  =1.0  MPaVin  (i.e.,  concrete). 

The  results  of  the  DEM  simulations  are  also  shown  in  Figure  5.5,  where  the 
hollow  data  points  represent  values  of  the  far  field  stress  that  cause  the  crack  to  propagate 
unstably  across  the  model.  As  mentioned  before,  using  Equation  4.3  for  this  problem 
gives  oF  =2.61  MPa  .  Applying  Richardson  extrapolation  to  the  results  of  these  three 
models  provides  an  estimate  of  oull  =2.595  MPa  as  the  value  of  the  far  field  stress  that 
the  DEM  models  will  converge  to,  which  is  within  0.5%  of  the  solution  in  Equation  4.3, 
which  could  be  in  error  by  as  much  as  0.5%.  Therefore,  the  proposed  criterion  for 
contact  failure  appears  very  promising. 

Applying  Equation  4.1  to  the  results  of  the  three  models,  a  convergence  rate 
p  =  1.8  (almost  quadratic)  was  obtained,  which  is  a  reasonable  value  and  is  encouraging. 
The  models  shown  in  Figure  5.5  were  also  developed  with  no  initial  crack,  and  all  three 
models  essentially  carried  a  maximum  stress  of  a  =  3.44  MPa  (independent  of  model 
discretization),  which  agrees  with  the  value  of  the  ultimate  tensile  strength  of  the  material 
crH/, .  Therefore,  the  proposed  method  provides  excellent  results  regardless  of  the  type  of 
failure  (fracture  toughness  or  ultimate  tensile  strength). 
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Figure  5.5  Convergence  for  crack  growth  in  a  DEM  model  of  a  cracked  plate  subjected 
to  pure  Mode  I  loading;  hollow  data  points  in  the  plot  represent  the  value  of  far  field 
stress  and  displacement  that  cause  unstable  crack  growth. 


6.  DEM  for  modeling  concrete  penetration 

One  of  the  objectives  of  our  study  was  to  use  DEM  for  modeling  impact  and 
penetration  events,  with  special  interest  on  concrete  members.  An  important  issue  is  the 
influence  of  strain  rate  effects  on  material  properties  such  as  the  tensile  strength  and 
elastic  modulus,  among  others.  Thus,  material  properties  determined  from  static  tests  may 
not  be  applicable  to  impact  situations  without  adjustments  to  account  for  the  rate  of 
loading.  In  a  study  of  fracture  of  concrete  under  uniaxial  impact  tensile  loading,  Zielinski 
(1982)  observed  that  the  elastic  stiffness,  stress  levels  at  failure,  and  deformation  capacity 
were  all  higher  when  compared  to  static  uniaxial  tensile  tests.  In  fact,  for  dynamic 
loading,  the  ultimate  compressive  strength  can  be  more  than  doubled  (Bischoff  and  Perry 
1991),  whereas  the  ultimate  uniaxial  tensile  strength  increases  by  multiples  of  5  to  7  at 
very  high  strain  rates  (Ross  et  al.  1996).  In  Tavarez  (2005),  we  review  some  of  the 
literature  on  strain  rate  effects  with  special  interest  in  concrete.  We  also  review  models 
that  have  been  proposed  by  Tang,  Malvern,  and  Jenkins  (1992),  Dilger,  Koch,  and 
Kowalczyk  (1984),  Mikkola  and  Sinisalo  (1982),  and  Soroushian,  Choi,  and  Alhamad 
(1986),  among  others. 

Tensile  Strength 

Far  less  dynamic  testing  has  been  carried  out  on  concrete  in  tension  than  in 
compression  (Fu  et  al.  1991),  and  therefore  there  is  insufficient  data  available  to  allow  the 
derivation  of  reliable  empirical  expressions.  Suaris  and  Shah  (1983)  suggested  that  rate 
sensitivity  is  considerably  more  pronounced  in  tension  and  flexure  than  in  compression. 
Tedesco  et  al.  (1989)  carried  out  split  Hopkinson  bar  tests  on  direct  tensile  specimens  and 
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observed  ratios  of  dynamic  to  static  strengths  of  up  to  3.6  at  strain  rates  of  approximately 
8  s'1.  Weerheim  and  Reinhardt  (1989)  carried  out  tensile  loading  of  notched  specimens 
using  a  split  Hopkinson  bar  and  found  a  ratio  of  dynamic  to  static  tensile  strength  of 
approximately  2  at  a  stress  rate  of  104  MPa/s,  rising  sharply  to  4.5  at  2xl05  MPa/s. 
Brara  et  al.  (2001)  conducted  experiments  on  concrete  specimens  in  tension  using  a  new 
experimental  technique  base  on  the  Hopkinson  bar  principle  and  the  spalling 
phenomenon,  and  concluded  that  for  strain  rates  of  approximately  120  s'1,  the  dynamic 
tensile  strength  exceeds  the  quasi-static  value  by  approximately  an  order  of  magnitude. 
John  et  al.  (1987)  proposed  a  model  to  predict  the  dynamic  tensile  strength  for  concrete 
using  three  material  parameters:  fracture  toughness  Klc,  critical  crack  tip  opening 
displacement  CTODc ,  and  Young’s  Modulus  E.  In  their  approach,  Klc,  and  E  are 
assumed  to  be  rate  independent  while  CTODc  is  assumed  to  decrease  exponentially  with 
the  logarithm  of  the  relative  strain  rate.  The  dynamic  tensile  strength  is  then  computed 
using  the  following  formula 


(O,-  1.4705  • 


E-{CTODc\, 
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where  the  dynamic  critical  crack  tip  opening  displacement  is  calculated  as 
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where  the  static  critical  crack  tip  opening  displacement  CTODc  =  0.0005  in.  and  the 
quasi-static  strain  rate  £0  =  30xl0'6  s'1. 


Elastic  Modulus 

The  limited  data  available  suggest  that  increases  in  Young’s  modulus  with  strain 
rate  are  far  less  pronounced  than  increases  in  compressive  strength.  Using  the  data 
available  from  the  literature,  Williams  (1994)  proposed  the  following  formula  for  the 
ratio  of  dynamic  to  static  Young’s  modulus  \EcdjEcs ) 

^  =  1.325 e  0  0234  £  a  6xl0'6  (6.4) 

E 
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However,  due  to  the  very  wide  scatter  of  the  data,  Williams  stated  that  this  relationship 
should  be  treated  with  caution.  The  following  section  will  discuss  the  actual  equations 
that  were  incorporated  into  our  DEM  method  to  treat  strain  rate  effects. 

Strain  Rate  Effects  Formulae  used  in  the  DEM  Method 

Although  more  work  needs  to  be  done  in  developing  models  for  rate  of  loading 
effects  for  concrete,  and  in  testing  applications  of  these  to  DEM  models,  we  have 
employed  the  following.  Klepaczco  (1990)  proposed  a  criterion  to  determine  the  dynamic 
tensile  strength  of  concrete  based  on  experimental  results  relating  fracture  stress  oF  and 
critical  loading  time  tc0 .  This  criterion  was  proposed  for  high  strain  rates  (higher  than  20 
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s  ')  and  has  the  advantage  of  being  based  on  some  physical  notion  taking  into  account  the 
thermally  activated  rate  process.  The  criterion  in  the  integral  form  is  expressed  as 
«(n 
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where  aull ,  tc0 ,  and  a{T)  are  three  material  constants  at  a  given  temperature  T.  The 
parameter  a  is  temperature  dependent  and  is  related  to  the  activation  energy  of  material 
separation.  The  parameter  tc0  is  the  longest  critical  time  at  which  oF(tcn)=  ouh  ,  where 
0U„  is  the  quasi-static  tensile  strength.  For  a  linearly  proportional  loading,  the  integral 
form  of  the  criterion  in  Equation  6.5  takes  the  following  form 


aF 


=  [(a  +  l)tc0Eau„ 


(6.6) 


where  E  is  the  Young’s  modulus  (under  static  conditions)  and  e  is  the  strain  rate. 
According  to  experimental  results  on  wet  concrete,  Brara  et.  al  (2001)  found  these 
parameters  to  be:  a  =  0.95  and  tc0  -  50  p.s.  Comparing  the  results  of  Equation  6.6  with 

data  obtained  from  the  literature,  it  was  found  that  the  formula  fits  very  well  to  the 
majority  of  experimental  data  and  can  also  be  used  for  smaller  strain  rate  values  (but 
greater  than  0.75  s'1)  and  still  obtain  reasonable  results.  As  such,  the  expression  in 
Equation  6.6  was  adopted  in  our  work  to  obtain  the  dynamic  tensile  strength  as  a  function 
of  strain  rate  for  strain  rates  greater  than  0.75  s'1.  For  smaller  strain  rates,  the  static  value 
of  tensile  strength  was  used.  In  our  DEM  code,  the  normal  strain  rate  is  computed  at 
each  time  step  using 

|v-n 
e  = J - 
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(6.7) 


where  v  is  the  relative  velocity  vector  between  two  DEM  particles,  n  is  the  normal  unit 
vector,  and  L0  is  the  initial  center- to-center  distance  between  bonded  DEM  elements. 


It  is  well  known  that  enhancement  of  elastic  modulus  due  to  strain  rate  effects  is 
less  significant  than  that  of  tensile  strength  and  compressive  strength  (Williams,  1994; 
Koh  et  al.  2001).  We  use  the  following  expressions  to  determine  the  dynamic  elastic 
modulus  as  a  function  of  strain  rate  in  compression  and  tension,  respectively  (CEB, 
1988): 

0.026 

(6.8) 
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where  the  subscripts  c,  t,  d,  and  s  denote  compressive,  tensile,  dynamic,  and  (quasi-) 
static  states,  respectively.  The  quasi-static  strain  rate  is  taken  as  e0  =  30xl0~6  s'1. 


It  has  also  been  observed  that  the  fracture  toughness  Klc  increases  as  the  rate  of 
applied  load  is  increased.  Mindess  et  al.  (1987)  conducted  experiments  on  single-edge 
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notched  concrete  beams  loaded  dynamically  in  3-point  bending,  using  an  instrumented 
drop  weight  impact  machine  in  order  to  study  the  effect  of  the  load  rate  on  the  fracture 
energy  and  fracture  toughness  of  the  material.  It  was  found  that  fracture  toughness 
values  and  fracture  energies  under  impact  loading  were  much  higher  (approximately  an 
order  of  magnitude)  than  those  obtained  in  static  tests. 

An  experimental  investigation  of  mode  I  fracture  of  concrete  was  conducted  by 
Lambert  and  Ross  (2000)  under  the  dynamic  loading  of  a  split  Hopkinson  pressure  bar. 
Fracture  specimens  in  the  form  of  notched-cavity  splitting  tension  cylinders  were 
subjected  to  stress  wave  loading  that  produced  strain  rates  up  to  10  s'1.  A  very  good 
linear  relationship  between  fracture  toughness  and  strain  rate  was  found  over  a  strain  rate 
range  of  1  s'1  to  8  s'1.  The  fracture  toughness  increased  by  a  factor  of  approximately  3 
when  the  strain  rate  increased  from  1  s'1  to  8  s'1.  However,  this  linear  trend  was  not 
suggested  as  a  global  trend  for  the  entire  strain  rate  spectrum.  It  is  important  to  mention 
that  the  analysis  of  fracture  toughness  data  under  high  rates  of  loading,  or  under  impact 
loading  is  complex.  For  instance,  the  dynamic  fracture  toughness  is  dependent  on  the 
crack  velocity,  which  may  vary  considerably  as  the  crack  propagates.  Therefore,  due  to 
the  lack  of  available  data  for  K,c  at  higher  strain  rates,  it  was  decided  instead  to  apply 
strain  rate  effects  to  the  material’s  fracture  energy  Gf,  since  the  failure  criterion  is 
primarily  based  on  this  quantity.  As  such,  it  was  also  assumed  that  the  material’s 
dynamic  fracture  energy  changes  in  the  same  fashion  as  the  dynamic  ultimate  tensile 
strength  as  a  function  of  strain  rate.  That  is,  we  used  the  following  expression  to 
determine  the  dynamic  fracture  energy  as  a  function  of  strain  rate 


(G/)„-HG/),  (6.10) 

where  the  factor  §  is  the  ratio  of  dynamic  to  static  tensile  strength  for  the  given  strain 
rate.  This  factor  is  directly  applied  to  the  computed  static  fracture  energy  (Equation  5.3) 
rather  than  to  K,c  so  that  strain  rate  effects  do  not  become  excessive  from  computing  the 
square  of  the  dynamic  fracture  toughness  in  order  to  determine  the  dynamic  fracture 
energy  using  the  relationship  in  Equation  5.3. 

DEM  Compression  Behavior  -  Material  Densification 

One  aspect  of  our  DEM  method  that  has  not  yet  been  described  is  the 
compression  behavior  of  the  medium.  Often,  the  possibility  of  compression  failures  in 
DEM  models  is  neglected,  with  the  main  emphasis  being  on  tensile  failures.  However,  in 
problems  involving  penetration,  extremely  high  compressive  stresses  are  produced,  and 
the  possibility  of  compressive  failure  should  be  included.  Experimental  results  have 
shown  that  concrete  targets  exhibit  what  is  referred  to  as  “medium  densification”  for  the 
material  ahead  of  the  projectile,  which  is  another  local  source  of  energy  absorption. 
There  is  likely  a  considerable  amount  of  energy  that  is  absorbed  due  to  this  phenomenon 
and  this  is  incorporated  to  some  extent  into  the  method  we  have  developed. 

To  illustrate  our  DEM  approach,  consider  the  bonded  elements  shown  in  Figure 
6.1(a).  The  cluster,  which  is  attached  to  a  boundary  that  could  be  a  fixed  support  or 
another  DEM  element,  consists  of  elements  of  initial  radius  r0 .  As  shown  in  Figure  6.1, 
the  initial  (equilibrium)  distance  between  the  elements  is  given  by  d0 .  Therefore,  contact 
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or  separation  is  determined  by  comparing  the  sum  of  the  radii  of  the  elements  with  the 
initial  distance  d0  (for  bonded  contacts). 

When  the  contact  is  in  compression,  the  normal-direction  spring  behaves  in  a 
linear  fashion  using  the  normal  spring  stiffness  given  in  Equation  3.1  for  plane  stress  and 
plane  strain.  However,  we  propose  that  when  the  contact  achieves  a  critical  compressive 
stress,  the  equilibrium  distance  d0  will  then  be  changed  to  d*Q ,  as  shown  in  Figure  6.1(b), 

and  the  contact  normal  force  will  behave  in  a  plastic  manner  in  order  to  account  for  non- 
recoverable  material  densification.  Therefore,  further  compression  (or  separation)  of  the 
contact  will  be  measured  relative  to  this  new  equilibrium  configuration  given  by  d*0 .  In 

essence,  the  DEM  elements  are  allowed  to  overlap.  Once  the  medium  has  reached  total 
compaction,  the  contact  will  behave  again  in  a  linear  fashion,  with  a  normal  contact 
stiffness  equal  to  its  initial  stiffness. 


Figure  6. 1  DEM  treatment  of  compressive  behavior  and  material  densification. 


Rather  than  allowing  DEM  elements  to  overlap,  the  densification  process  described  here 
could  also  be  implemented  by  allowing  the  element  radius  to  reduce  from  r0  to  r\  as 
suggested  in  Figure  6.1(b),  and  increasing  the  element  density  in  order  to  conserve  the 
element’s  mass.  However,  this  latter  approach  was  not  used  in  our  study  because 
reduction  of  the  element  radius  will  most  likely  affect  other  contacts  associated  with  this 
element  (i.e.,  contact  separation  will  occur,  as  shown  in  the  figure).  In  the  approach  we 
use  wherein  some  DEM  elements  are  allowed  to  overlap,  other  contacts  associated  with 
the  element  (which  may  not  be  at  a  critical  compressive  stage)  will  be  unchanged,  and 
material  densification  will  automatically  be  obtained  because  of  the  overlap  between 
appropriate  DEM  elements.  This  technique  is  also,  in  some  aspects,  analogous  to 
establishing  a  nonlinear  equation  of  state  (EOS)  for  the  medium. 

In  order  to  determine  the  compressive  critical  contact  stress,  the  results  obtained 
from  the  DEM  unit  cell  shown  in  Figure  5.3  were  used.  The  critical  normal  force  given 
in  Equation  5.5  is  then  applied,  using  the  value  of  maximum  compressive  strength  fc 
instead  of  the  maximum  tensile  strength  oull .  The  critical  compressive  contact  force 
(f„cr“  )  is  then  given  by 


(6.11) 
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Since  the  material’s  compressive  strength  is  now  included,  this  parameter  will  also 
be  dependent  on  the  strain  rate  i  .  The  dynamic  compressive  strength  will  be  determined 
using  (CEM,  1988) 
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(6.13) 


where 


log  y  =  6.15a  -  0.492 


(6.14) 


a -1/(5  + 3/; /4) 


(6.15) 


and  fc  is  the  static  compressive  strength  in  MPa  at  a  quasi-static  strain  rate 
£0  =  30xl(T6  s'1.  Figure  6.2  shows  the  ratio  of  dynamic  to  static  compressive  strength 
of  concrete  as  a  function  of  strain  rate.  The  strain  rate  has  a  smaller  effect  on 
compressive  strength  than  on  the  tensile  strength  of  the  material. 


strain  rats,  s'1 

Figure  6.2  Strain  rate  effect  on  ultimate  compressive  strength. 

It  is  important  to  mention  that  due  to  the  lack  of  available  data  for  validation,  it 
was  not  clear  how  to  choose  or  determine  the  parameters  that  govern  the  compaction 
phase.  These  parameters  are  the  compaction  stiffness  and  the  strain  at  full  compaction 
and  they  were  chosen  somewhat  arbitrarily.  As  mentioned  before,  is  was  assumed  that  in 
the  compaction  phase  the  contact  behaves  plastically  and  it  will  continue  to  do  so  until 
the  distance  between  the  element  centers  (given  by  d 0* )  reaches  the  DEM  element  radius 
R .  When  this  happens,  it  is  assumed  that  the  material  has  reached  full  compaction  and 
therefore,  it  will  continue  to  have  a  normal  stiffness  equal  to  the  initial  stiffness  Kn . 
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Additional  work  is  required  to  validate  the  parameters  chosen  for  this  compaction  phase 
and  to  gain  more  confidence  in  their  use.  Nevertheless,  as  will  be  shown  later,  it  is  an 
essential  phenomenon  to  incorporate  in  a  high  speed  penetration  event,  especially  for 
concrete  members. 


7.  DEM  simulation  of  concrete  penetration 

This  work  developed  a  new  approach  for  using  the  DEM  method  to  assess  local 
damage  to  concrete  structures  subjected  to  extreme  loading.  Currently,  damage  to 
reinforced  concrete  structures  is  evaluated  mainly  with  empirical  formulas  derived  from 
various  types  of  impact  tests.  DEM  simulation  results  are  presented  here  with  a 
comparison  with  available  formulas  for  predicting  penetration,  scabbing  thickness,  and 
perforation  thickness  in  plain  concrete  targets.  A  limitation  to  our  numerical  DEM 
implementation  is  that  real  life  penetration  events  are  three  dimensional,  whereas  our 
DEM  model  is  presently  two  dimensional.  Tavarez  (2005)  describes  a  calibration  process 
so  that  qualitative  comparisons  can  be  made  between  our  two  dimensional  simulations 
and  the  available  empirical  formulas.  We  note  that  for  two  dimensional  simulations,  no 
calibration  process  is  needed,  but  unfortunately,  data  for  structure  behavior  under 
extreme  loading  conditions  is  not  as  well  developed.  For  more  accurate  applications  to 
projectile  penetration,  a  three  dimensional  implementation  of  our  results  is  needed,  and 
this  may  be  considered  in  future  work. 

Figure  7.1  shows  the  DEM  model  used  for  simulating  penetration  events.  The 
concrete  member  consists  of  12000  DEM  elements  and  is  cantilever-supported  at  both 
ends.  The  model  has  a  length  of  96  in.  with  a  depth  of  30  in.,  and  is  assumed  to  be  in 
plane  strain  conditions.  Therefore,  all  quantities  (density,  contact  stiffness,  stresses,  etc.) 
are  computed  per  unit  thickness.  The  projectile  is  steel  and  consists  of  351  DEM 
elements  and  was  constructed  using  a  close-packed  assembly.  The  projectile  has  a  width 
(diameter)  of  3.5  in.  and  its  density  was  determined  so  that  its  total  weight  is  15  lbs.  The 
following  subsections  provide  more  details  about  the  properties  used  for  both  materials. 


|- -  96  In  - H 

Figure  7.1  DEM  discretization  for  modeling  penetration  in  a  concrete  member. 

Concrete  Target 

Table  7. 1  shows  the  DEM  model  details  and  material  properties  for  the  concrete 
target  used  in  all  simulations.  All  material  properties  reported  in  this  table  are  quasi-static 
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values.  Such  values  are  then  used  in  the  code  to  compute  the  respective  dynamic 
properties  as  a  function  of  strain  rate. 


Table  7. 1  Concrete  properties  used  in  DEM  model. 


Density 

(p) 

DEM  element 
radius 
(K) 

Elastic 

modulus 

(E) 

Poisson's 

ratio 

M 

Tensile 

strength 

(O’ oft) 

Compressive 

strength 

(f'c) 

Fracture 
toughness 
(K, c) 

2.17x10"'  lb-s2/in 

0.25  in. 

3,122  ksi 

0.2 

300  psi 

3000  psi 

910  psi-s1/2 

Steel  Projectile 

Table  7.2  shows  the  DEM  model  details  and  material  properties  used  for  the  steel 
projectile  used  in  the  analysis.  It  is  important  to  mention  that  even  though  the  projectile 
is  deformable,  it  does  not  have  failure  capabilities  incorporated  for  it.  This  was  decided 
due  to  the  relative  coarse  mesh  (resulting  in  a  somewhat  rough  surface)  and  assembly 
used  for  the  projectile  model.  Moreover,  experimental  data  show  that  most  projectile 
damage  phenomena  is  due  to  thermodynamic  issues  and  therefore,  since  this  aspect  of 
deformation  was  not  incorporated  into  the  model,  it  was  decided  not  to  include  erosion  or 
failure  of  the  projectile. 


Table  7.2  Projectile  properties  used  in  DEM  model. 


Effective 

Density 

(P) 

DEM  element 
radius 
(R) 

Elastic 

modulus 

(£) 

Poisson’s 

ratio 

(y) 

Projectile 

diameter 

(D) 

Projectile 

weight 

( w ) 

2.43x1  O'3  lb-s2/in 

0.15544  in 

29,000  ksi 

0.3 

3.5  in 

151b 

DEM  Analysis  Procedure 

The  velocity  of  the  projectile  was  varied  in  the  numerical  simulations,  and 
depending  on  the  impact  velocity,  different  behaviors  were  observed.  However,  in  order 
to  compare  the  numerical  results  produced  by  the  2-D  DEM  code  used  in  this  work  with 
available  empirical  formulas  for  penetration,  scabbing,  and  perforation  of  concrete 
targets,  which  are  based  on  3-D  experimental  results,  the  following  approach  was  taken. 
After  creating  the  DEM  model  shown  in  Figure  7.1,  it  was  necessary  to  calibrate  the 
model  to  determine  the  effective  tensile  strength  and  fracture  toughness  of  the  medium. 
According  to  the  available  empirical  formulas  and  the  given  material  properties,  for  the 
model  considered  here,  the  penetration  depth  should  be  approximately  5  inches  for  a 
projectile  velocity  of  500  ft/s,  without  producing  scabbing  in  the  back  face  of  the  target. 
As  such,  the  initial  velocity  of  the  projectile  in  the  DEM  model  was  set  to  V0  =  500  ft/s, 
and  the  calibration  parameter  was  varied  until  a  penetration  depth  of  approximately  5 
inches  was  obtained.  Once  this  calibration  was  completed,  subsequent  simulations  were 
carried  out  with  no  recalibrations.  In  the  subsequent  simulations: 

•  The  velocity  of  the  projectile  was  gradually  increased  in  order  to  produce  scabbing, 

and  then  penetration  in  the  concrete  target. 

•  The  scabbing  and  penetration  velocities  obtained  from  the  numerical  simulations 

were  then  compared  to  the  results  obtained  using  the  available  empirical  formulae. 
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According  to  the  DEM  model  parameters,  the  required  time  step  for  the  simulations  was 
approximately  At  =  8.9 xlO-8  seconds.  Due  to  the  number  of  elements  in  the  model,  an 
average  running  time  of  approximately  15  hours  on  a  Unix  workstation  was  needed  for  a 
simulation  time  of  approximately  10  ms,  which  was  the  simulation  time  required  for  a 
complete  penetration  event. 


DEM  Simulation  Results 

Figure  7.2(a)  shows  the  projectile  at  its  maximum  depth  of  penetration,  which  is 
approximately  5.05  in.,  and  Figure  7.2(b)  shows  the  vertical  position  of  the  projectile  tip 
as  a  function  of  time.  The  hollow  data  point  indicates  when  the  projectile  impacts  the  top 
of  the  beam. 


lime,  ms 


Figure  7.2  Projectile  maximum  penetration  and  tip  vertical  position  as  a  function  of 
time. 


The  figure  also  shows  the  depth  of  the  beam,  which  is  shown  in  the  scale  of  the  figure 
(i.e.,  30  in.).  As  seen  in  Figure  7.2(a),  there  were  flexural  cracks  resulting  in  global 
deformation  in  the  target,  which  makes  it  difficult  to  determine  the  actual  penetration 
depth  relative  to  the  target.  Nevertheless,  the  results  obtained  in  this  simulation  agree 
very  well  with  the  NDRC  prediction. 

Figure  7.3  shows  the  results  of  the  different  available  empirical  formulae  to 
determine  the  scabbing  thickness  hs  of  a  concrete  target  subjected  to  projectile 
penetration.  As  defined  earlier,  the  scabbing  thickness  is  the  minimum  thickness  required 
to  prevent  scabbing  in  a  concrete  target. 

Since  it  is  a  lengthy  process  to  construct  DEM  models  of  different  thickness,  it 
was  instead  decided  to  determine  the  minimum  projectile  velocity  required  to  induce 
scabbing  in  the  model.  In  Figure  7.3,  we  see  that  for  the  30  in.  target  thickness  of  our 
DEM  model,  the  predicted  values  for  scabbing  thickness  vary  from  750  ft/s  to  1250  ft/s. 
In  order  to  determine  the  “scabbing  velocity”  in  our  DEM  model,  the  projectile  velocity 
was  increased  in  intervals  of  50  ft/s  until  scabbing  was  produced  in  the  model. 
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Figure  7.3  Calculated  scabbing  thickness  (using  data  listed  in  Tables  7. 1  and  7.2). 

Figure  7.4  shows  scabbing  in  the  DEM  model  for  a  projectile  impact  velocity  of 
V0  =  850  ft/s.  There  was  no  scabbing  at  a  velocity  of  800  ft/s  (only  a  large  flexural  crack 
that  leaded  to  ultimate  fracture  of  the  beam),  and  therefore  it  was  determined  that  the 
“scabbing  velocity”  for  the  model  is  850  ft/s  which  falls  within  the  predicted  values  of 
the  available  empirical  formulas.  The  snapshot  at  Stage  4  in  Figure  8.8  shows  the 
maximum  penetration  in  the  model. 


Figure  7.4  Scabbing  in  concrete  target. 

The  scabbing  velocity  for  the  model  agrees  very  well  with  the  prediction  of  the  BRL 
formula,  which  predicts  a  scabbing  velocity  of  V0  =  760  ft/s.  We  note  that  this  equation 
is  the  only  one  among  the  many  empirical  formulas  that  does  not  rely  on  the  penetration 
depth  Xj  to  determine  the  scabbing  thickness.  According  to  the  BRL  equations,  the 
scabbing  thickness  is  determined  from  the  perforation  thickness.  It  is  important  to 
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mention  that  the  addition  of  the  nonlinear  compression  behavior  decreases  the  degree  of 
brittleness  in  the  model  and  therefore  scabbing  is  produced  at  a  higher  velocity  than 
might  actually  occur.  As  seen  in  Figure  7.4  (stage  4),  considerable  penetration  occurred 
in  the  model  before  any  scabbing  was  present,  suggesting  that  when  it  did  occur,  there 
were  other  factors  that  contributed  to  it  such  as  overall  deformation  and  the  development 
of  flexural  and  shear  cracks. 

Figure  7.5  shows  the  results  of  the  different  available  empirical  formulae  to 
determine  the  perforation  thickness  hp  of  a  concrete  target,  where  the  perforation 

thickness  is  the  minimum  thickness  required  to  prevent  perforation  in  a  concrete  target. 
Analogous  to  the  procedure  done  for  scabbing,  it  was  decided  to  determine  the  minimum 
velocity  required  to  completely  perforate  the  target.  This  velocity  is  also  known  as  the 
ballistic  limit.  In  Figure  7.5,  we  can  see  that  for  a  target  thickness  of  30  in.,  the  predicted 
values  for  perforation  thickness  vary  from  approximately  1000  ft/s  to  1500  ft/s.  Figure 
7.6  shows  a  DEM  model  where  full  perforation  was  obtained.  As  in  the  procedure  done 
to  determine  the  scabbing  velocity,  the  projectile  velocity  was  increased  in  intervals  of  50 
ft/s  until  full  perforation  was  achieved  in  the  model.  This  was  obtained  for  an  impact 
velocity  of  1150  ft/s,  which  again  falls  within  the  range  of  predictions  from  the  available 
empirical  formulas. 

The  perforation  velocity  obtained  from  the  DEM  simulations  agrees  particularly 
well  with  the  prediction  from  the  BRL  equation  (shown  in  Figure  7.5),  which  directly 
predicts  the  perforation  thickness  without  using  the  penetration  depth  xf  or  the  scabbing 

thickness  hp.  Although  the  results  obtained  from  the  simulations  presented  in  this 

chapter  are  very  encouraging,  additional  work  is  required  so  that  more  confidence  is 
gained  in  the  use  of  the  proposed  approach  to  predict  local  damage  in  concrete  targets 
due  to  projectile  impact. 


projectile  velocity,  ft/s 

Figure  7.5  Calculated  perforation  thickness  (using  data  listed  in  Tables  7.1  and  7.2). 
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